Age Estimates
#......................
# read in data
#......................
modout <- readRDS("results/ModFits/GBR_age_rung50_burn10000_smpl10000.RDS")
Log Plot
modout$mcmcout$output %>%
dplyr::filter(rung == "rung1") %>%
dplyr::filter(stage == "sampling") %>%
dplyr::select(c("chain", "iteration", "loglikelihood", "logprior")) %>%
tidyr::gather(., key = "like", value = "val", 3:4) %>%
ggplot() +
geom_line(aes(x = iteration, y = val, color = like), size = 0.1, alpha = 0.8) +
facet_wrap(chain~like, scales = "free_y") +
theme_bw()
Quick MCMC Diagnostics
quickout <- quick_mc_diagnostics(modout)
quickout
IFRs & Incidence Curve
agekey <- modout$inputs$IFRmodel$IFRdictkey %>%
dplyr::rename(param = Strata)
#......................
# get ifrs
#......................
ifrs <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, whichrung = paste0("rung", 1),
what = "IFRparams", by_chain = FALSE)
ifrs <- dplyr::left_join(agekey, ifrs, by = "param") %>%
dplyr::mutate(param = factor(param, levels = paste0("ma", 1:5))) %>%
dplyr::arrange()
#......................
# get crude ifrs
#......................
dscdat <- readRDS("data/derived/descriptive_results_datamap.RDS")
GBRcrudedat <- dscdat %>%
dplyr::filter(study_id == "GBR2" & breakdown == "ageband") %>%
dplyr::select(c("study_id", "plotdat")) %>%
tidyr::unnest(cols = "plotdat")
ageplotdat <- GBRcrudedat %>%
dplyr::filter(seromidpt == obsday) %>%
dplyr::mutate(infxns = popn * seroprev,
crudeIFR = cumdeaths/infxns) %>%
dplyr::select(c("ageband", "crudeIFR"))
#......................
# get incidence curve
#......................
infxncurve <- COVIDCurve::draw_posterior_infxn_cubic_splines(IFRmodel_inf = modout,
dwnsmpl = 1e3,
by_chain = FALSE,
by_strata = TRUE)
#......................
# make ifr and incidence plots
#......................
plot1 <- ggplot() +
geom_pointrange(data = ifrs, aes(x = ageband, ymin = LCI, ymax = UCI, y = median),
color = "#969696", size = 1.2) +
geom_point(data = ageplotdat, aes(x = ageband, y = crudeIFR), color = "#000000", size = 2) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"),
legend.position = "right") +
xlab("") + ylab("Median (95% CIs)") +
labs(caption = "Grey and black points represent model fits and crude IFR estimates, respectively")
plot2 <- infxncurve$plotObj
# main plot
main_plotObj <- cowplot::plot_grid(plot1, plot2, ncol = 1, nrow = 2)
# quick out
jpgsnapshot("figures/GBR2_ifr_plot_age.jpg", plot = main_plotObj)
plot(main_plotObj)
ifrs %>%
dplyr::mutate_if(is.numeric, round, 2) %>%
pretty_DT_tab(.)
Overall IFR
# population denom
demog <- modout$inputs$IFRmodel$demog %>%
dplyr::rename(param = Strata)
demog <- demog %>%
dplyr::left_join(., agekey)
gIFR <- sum(ifrs$median * demog$popN/sum(demog$popN))
tibble::tibble("Overall IFR" = sprintf("%1.2f%%", 100*gIFR)) %>%
knitr::kable(.)
Seroprevalence Observed vs. Inferred True Prev
seropnts <- COVIDCurve::draw_posterior_sero_curves(IFRmodel_inf = modout,
dwnsmpl = 1e3,
by_chain = F)
serocurvedat <- seropnts %>%
dplyr::select(c("sim", "ObsDay", dplyr::starts_with("RG_"), dplyr::starts_with("inf_"))) %>%
tidyr::gather(., key = "seroprev_strata_lvl", value = "serocount", 3:ncol(.)) %>%
dplyr::mutate(seroprevlvl = ifelse(stringr::str_detect(seroprev_strata_lvl, "RG_"), "RG", "IF"),
param = stringr::str_extract(seroprev_strata_lvl, "ma[0-9]+")) %>%
dplyr::group_by_at(c("sim", "ObsDay", "seroprevlvl", "param")) %>%
dplyr::summarise(serocount = sum(serocount, na.rm = T) ) %>%
dplyr::left_join(., agekey, by = "param") %>%
dplyr::left_join(., demog, by = "ageband") %>%
dplyr::mutate(seroprev = serocount/popN)
# crude data
SeroPrevObs <- GBRcrudedat %>%
dplyr::select(c("obsdaymin", "obsdaymax", "ageband", "seroprev")) %>%
dplyr::filter(!duplicated(.)) %>%
dplyr::left_join(., demog, by = "ageband") %>%
dplyr::mutate(obsmidday = (obsdaymin + obsdaymax)/2)
serocurvedat %>%
ggplot() +
geom_line(aes(x = ObsDay, y = seroprev, color = seroprevlvl), alpha = 0.5) +
geom_rect(data = SeroPrevObs, aes(xmin = obsdaymin, xmax = obsdaymax, ymin = -Inf, ymax = Inf),
color = "#d9d9d9", fill = "#d9d9d9", alpha = 0.8) +
geom_point(data = SeroPrevObs, aes(x = obsmidday, y = seroprev, group = ageband),
color = "#000000", size = 1.2) +
facet_wrap(.~ageband) +
scale_color_manual("Seroprev. \n Adjustment", values = c("#FFD301", "#246BCF"),
labels = c("Inferred 'Truth'", "Inferred 'Observed' - \n Rogan-Gladen Corrected")) +
labs(caption = "Grey box is observed seroprevalence across study period") +
xyaxis_plot_theme +
theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"))
Posterior Predictive Check
#......................
# get posterior draws
#......................
postdat <- COVIDCurve::posterior_check_infxns_to_death(IFRmodel_inf = modout,
dwnsmpl = 1e3,
by_chain = FALSE)
postdat_long <- postdat %>%
dplyr::select(c("sim", "time", dplyr::starts_with("deaths"))) %>%
tidyr::gather(., key = "param", value = "deaths", 3:ncol(.)) %>%
dplyr::mutate(param = gsub("deaths_", "", param)) %>%
dplyr::left_join(., y = agekey, by = "param")
#......................
# get deaths posterior pred check
#......................
datclean <- modout$inputs$IFRmodel$data$obs_deaths %>%
dplyr::filter(Deaths != -1) %>%
dplyr::left_join(., modout$inputs$IFRmodel$IFRdictkey, by = "Strata")
#......................
# make plot
#......................
ppc_plot_age <- ggplot() +
geom_line(data = postdat_long, aes(x= time, y = deaths, group = ageband),
size = 1.2, color = "#bdbdbd") +
geom_line(data = datclean, aes(x=ObsDay, y = Deaths, group = ageband),
color = "#3182bd") +
facet_wrap(.~ageband) +
theme_bw() +
ggtitle("Posterior Predictive Check", subtitle = "Grey Lines are Draws from Posterior, Blue Line is Observed Data (from ECDC)") +
theme(plot.title = element_text(hjust = 0.5, vjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, vjust = 0.5))
jpgsnapshot("figures/GBR2_ppc_plot_age.jpg", plot = ppc_plot_age)
plot(ppc_plot_age)
Regional Estimates
#......................
# read in data
#......................
modout <- readRDS("results/ModFits/GBR_rgn_rung50_burn10000_smpl10000.RDS")
Log Plot
modout$mcmcout$output %>%
dplyr::filter(rung == "rung1") %>%
dplyr::filter(stage == "sampling") %>%
dplyr::select(c("chain", "iteration", "loglikelihood", "logprior")) %>%
tidyr::gather(., key = "like", value = "val", 3:4) %>%
ggplot() +
geom_line(aes(x = iteration, y = val, color = like), size = 0.1, alpha = 0.8) +
facet_wrap(chain~like, scales = "free_y") +
theme_bw()
Quick MCMC Diagnostics
quickout <- quick_mc_diagnostics(modout)
quickout
IFRs & Incidence Curve
rgnkey <- modout$inputs$IFRmodel$IFRdictkey %>%
dplyr::rename(param = Strata)
#......................
# get ifrs
#......................
ifrs <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, whichrung = paste0("rung", 1),
what = "IFRparams", by_chain = FALSE)
ifrs <- dplyr::left_join(rgnkey, ifrs, by = "param") %>%
dplyr::mutate(param = factor(param, levels = paste0("ma", 1:3))) %>%
dplyr::arrange()
#......................
# get crude ifrs
#......................
dscdat <- readRDS("data/derived/descriptive_results_datamap.RDS")
GBRcrudedat <- dscdat %>%
dplyr::filter(study_id == "GBR2" & breakdown == "region") %>%
dplyr::select(c("study_id", "plotdat")) %>%
tidyr::unnest(cols = "plotdat")
rgnplotdat <- GBRcrudedat %>%
dplyr::filter(seromidpt == obsday) %>%
dplyr::mutate(infxns = popn * seroprev,
crudeIFR = cumdeaths/infxns) %>%
dplyr::select(c("region", "crudeIFR"))
#......................
# get incidence curve
#......................
infxncurve <- COVIDCurve::draw_posterior_infxn_cubic_splines(IFRmodel_inf = modout,
dwnsmpl = 1e3,
by_chain = FALSE,
by_strata = TRUE)
#......................
# make ifr and incidence plots
#......................
plot1 <- ggplot() +
geom_pointrange(data = ifrs, aes(x = region, ymin = LCI, ymax = UCI, y = median),
color = "#969696", size = 1.2) +
geom_point(data = rgnplotdat, aes(x = region, y = crudeIFR), color = "#000000", size = 2) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"),
legend.position = "right") +
xlab("") + ylab("Median (95% CIs)") +
labs(caption = "Grey and black points represent model fits and crude IFR estimates, respectively")
plot2 <- infxncurve$plotObj
# main plot
main_plotObj <- cowplot::plot_grid(plot1, plot2, ncol = 1, nrow = 2)
# quick out
jpgsnapshot("figures/GBR2_ifr_plot_rgn.jpg", plot = main_plotObj)
plot(main_plotObj)
ifrs %>%
dplyr::mutate_if(is.numeric, round, 2) %>%
pretty_DT_tab(.)
Overall IFR
# population denom
demog <- modout$inputs$IFRmodel$demog %>%
dplyr::rename(param = Strata)
demog <- demog %>%
dplyr::left_join(., rgnkey)
gIFR <- sum(ifrs$median * demog$popN/sum(demog$popN))
tibble::tibble("Overall IFR" = sprintf("%1.2f%%", 100*gIFR)) %>%
knitr::kable(.)
Seroprevalence Observed vs. Inferred True Prev
seropnts <- COVIDCurve::draw_posterior_sero_curves(IFRmodel_inf = modout,
dwnsmpl = 1e3,
by_chain = F)
serocurvedat <- seropnts %>%
dplyr::select(c("sim", "ObsDay", dplyr::starts_with("RG_"), dplyr::starts_with("inf_"))) %>%
tidyr::gather(., key = "seroprev_strata_lvl", value = "serocount", 3:ncol(.)) %>%
dplyr::mutate(seroprevlvl = ifelse(stringr::str_detect(seroprev_strata_lvl, "RG_"), "RG", "IF"),
param = stringr::str_extract(seroprev_strata_lvl, "ma[0-9]+")) %>%
dplyr::group_by_at(c("sim", "ObsDay", "seroprevlvl", "param")) %>%
dplyr::summarise(serocount = sum(serocount, na.rm = T) ) %>%
dplyr::left_join(., rgnkey, by = "param") %>%
dplyr::left_join(., demog, by = "region") %>%
dplyr::mutate(seroprev = serocount/popN)
# crude data
SeroPrevObs <- GBRcrudedat %>%
dplyr::select(c("obsdaymin", "obsdaymax", "region", "seroprev")) %>%
dplyr::filter(!duplicated(.)) %>%
dplyr::left_join(., demog, by = "region") %>%
dplyr::mutate(obsmidday = (obsdaymin + obsdaymax)/2)
serocurvedat %>%
ggplot() +
geom_line(aes(x = ObsDay, y = seroprev, color = seroprevlvl), alpha = 0.5) +
geom_rect(data = SeroPrevObs, aes(xmin = obsdaymin, xmax = obsdaymax, ymin = -Inf, ymax = Inf),
color = "#d9d9d9", fill = "#d9d9d9", alpha = 0.8) +
geom_point(data = SeroPrevObs, aes(x = obsmidday, y = seroprev, group = region),
color = "#000000", size = 1.2) +
facet_wrap(.~region) +
scale_color_manual("Seroprev. \n Adjustment", values = c("#FFD301", "#246BCF"),
labels = c("Inferred 'Truth'", "Inferred 'Observed' - \n Rogan-Gladen Corrected")) +
labs(caption = "Grey box is observed seroprevalence across study period") +
xyaxis_plot_theme +
theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"))
Posterior Predictive Check
#......................
# get posterior draws
#......................
postdat <- COVIDCurve::posterior_check_infxns_to_death(IFRmodel_inf = modout,
dwnsmpl = 1e3,
by_chain = FALSE)
postdat_long <- postdat %>%
dplyr::select(c("sim", "time", dplyr::starts_with("deaths"))) %>%
tidyr::gather(., key = "param", value = "deaths", 3:ncol(.)) %>%
dplyr::mutate(param = gsub("deaths_", "", param)) %>%
dplyr::left_join(., y = rgnkey, by = "param")
#......................
# get deaths posterior pred check
#......................
datclean <- modout$inputs$IFRmodel$data$obs_deaths %>%
dplyr::filter(Deaths != -1) %>%
dplyr::left_join(., modout$inputs$IFRmodel$IFRdictkey, by = "Strata")
#......................
# make plot
#......................
ppc_plot_rgn <- ggplot() +
geom_line(data = postdat_long, aes(x= time, y = deaths, group = region),
size = 1.2, color = "#bdbdbd") +
geom_line(data = datclean, aes(x=ObsDay, y = Deaths, group = region),
color = "#3182bd") +
facet_wrap(.~region) +
theme_bw() +
ggtitle("Posterior Predictive Check", subtitle = "Grey Lines are Draws from Posterior, Blue Line is Observed Data (from ECDC)") +
theme(plot.title = element_text(hjust = 0.5, vjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, vjust = 0.5))
jpgsnapshot("figures/GBR2_ppc_plot_rgn.jpg", plot = ppc_plot_rgn)
plot(ppc_plot_rgn)
